The present invention is generally directed to a system and method for recovering wave front phase information and using the recovered information. More particularly, the invention is directed to a system and method for determining the phase information associated with a wave front from measured intensity information.
Huygens, Kirchhoff, Sommerfield and Rayleigh originated and contributed most to the currently accepted theory of diffraction, which forms the theoretical foundation for the present invention. Basically, the theory postulates that given a known wave front filling a planar window in an otherwise dark planar screen, the wave front at any point beyond the screen is calculable. Variations of this theory are used to compute the Fraunhofer far-field antenna pattern assuming a known field distribution at the antenna for electromagnetic wave fronts in the microwave range. An ordinary light camera, assuming a sufficiently coherent and pseudo monochromatic light wave, contains the Back Focal Plane (BFP) of the camera""s objective lens as the illuminating window of the diffraction theory, and the imaging plane as the plane at which the image could be calculated. Of course, in the case of the camera, photographic film or electronic sensing devices are placed in the image plane, recording the intensity of the wave and no calculations need be made. However, it will be appreciated that at each point in a wave front there is more than just the intensity of the wavexe2x80x94there is the phase of the wave which may contain as much as eighty percent of the information about the object which is being imaged. To appreciate this fact more fully, it is only necessary to recall the potential of the conventional holograms to image objects in three dimensions. In particular, using phase information about a coherent wave front, holography creates three-dimensional images such that obscured objects may become visible if the observer steps to the side. So, the problem addressed using this invention can be stated as follows: given that a wave front is a complex function characterized by both amplitude (related in a straightforward manner to intensity) and phase at each point, how can the phase be captured using only intensity measurements.
To appreciate the complexity of the problem, consider the following observation: at a given instant of time, the phase of a wave length is about 6.28 radians. For red light, the length over which that phase is generated is about 0.6 microns. Considering that light travels at approximately 300,000,000 meters per second, that means that the frequency of such a wave passing a point in space is about 3.1*1015 radians/second. No device exists that has that kind of response time. For the ordinary light camera, the two planes of interest relating to diffraction theory are the BFP of the lens and the image plane. They have been shown to be conjugate planes in the sense that the wave front in the image plane is essentially the Fourier Transform of the illuminating wave in the BFP.
In a coherent monochromatic imaging system the problem of extracting phase information from a detection medium which records only intensity information remains a problem without a consistent solution. Several experimental methods have been proposed for determining the phase function across a wave front. One such method disclosed in Gabor, D. xe2x80x9cA New Microscope Principle,xe2x80x9d Nature 161, 777 (1948) involves the addition of a reference wave to the wave of interest in the recording plane. The resulting hologram records a series of intensity fringes, on a photographic plate, which contain enough information to reconstruct the complete wave function of interest. However, in most practical applications this method is cumbersome and impractical to employ.
Other methods, which do not employ reference waves, have been proposed for inferring the complete wave function from intensity recordings. See, e.g., Erickson, H. and Klug, A. xe2x80x9cThe Fourier Transform of an Electron Micrograph: Effects of Defocusing and Aberrations, and Implications for the use of Underfocus Contrast Enhancementsxe2x80x9d, Berichte der Bunsen Gesellschaft, Bd. 74, Nr. 11, 1129-1137 (1970). For the most part, these methods involve linear approximation and thus are only valid for small phase and/or amplitude deviations across the wave front of interest. In general, these methods also suffer from the drawback of requiring intensive computational resources.
A further method proposed that intensity recordings of wave fronts can be made conveniently in both the imaging and diffraction planes. Gerchberg, R. and Saxton, W. xe2x80x9cPhase Determination from Image and Diffraction Plane Pictures in the Electron Microscope,xe2x80x9d Optik, Vol. 34, No. 3, pp. 275-284 (1971). The method uses sets of quadratic equations that define the wave function across the wave in terms of its intensity in the image and diffraction planes. This method of analysis is not limited by the above-described deficiency of being valid for small phase or amplitude deviations, but again, in general it requires a large amount of computational resources.
In 1971 the present inventor co-authored a paper describing a computational method for determining the complete wave function (amplitudes and phases) from intensity recordings in the imaging and diffraction planes. See, xe2x80x9cA Practical Algorithm for the Determination of Phase from Image and Diffraction Plane Pictures,xe2x80x9d Cavendish Laboratory, Cambridge, England, Optik, Vol. 35, No. 2, (1972) pp. 237-246, which is incorporated herein by reference for background. The method depends on there being a Fourier Transform relation between the complex wave functions in these two planes. This method has proven to have useful applications in electron microscopy, ordinary light photography and crystallography where only an x-ray diffraction pattern may be measured.
The so-called Gerchberg-Saxton solution is depicted in a block diagram form in FIG. 1. The input data to the algorithm are the square roots of the physically sampled wave function intensities in the image 100 and diffraction 110 planes. Although instruments can only physically measure intensities, the amplitudes of the complex wave functions are directly proportional to the square roots of the measured intensities. A random number generator is used to generate an array of random numbers 120 between xcfx80 and xe2x88x92xcfx80, which serve as the initial estimates of the phases corresponding to the sampled imaged amplitudes. If a better phase estimate is in hand a priori, that may be used instead. In step 130 of the algorithm, the estimated phases 120 (represented as unit amplitude xe2x80x9cphasorsxe2x80x9d) are then multiplied by the corresponding sampled image amplitudes from the image plane, and the Discrete Fourier Transform of the synthesized complex discrete function is accomplished in step 140 by means of the Fast Fourier Transform (FFT) algorithm. The phases of the discrete complex function resulting from this transformation are retained as unit amplitude xe2x80x9cphasorsxe2x80x9d (step 150), which are then multiplied by the true corresponding sampled diffraction plane amplitudes in step 160. This discrete complex function (an estimate of the complex diffraction plane wave) is then inverse Fast Fourier transformed in step 170. Again the phases of the discrete complex function generated are retained as unit amplitude xe2x80x9cphasorsxe2x80x9d (step 180), which are then multiplied by the corresponding measured image amplitudes to form the new estimate of the complex wave function in the image plane 130. The sequence of steps 130-180 is then repeated until the computed amplitudes of the wave forms match the measured amplitudes sufficiently closely. This can be measured by using a fraction whose numerator is the sum over all sample points in either plane of the difference between the measured and computed amplitudes of the complex discrete wave function squared and whose denominator is the sum over all points in the plane of the measured amplitudes squared. When this fraction is less than 0.01 the function is usually well in hand. This fraction is often described as the sum of the squared error (SSE) divided by the measured energy of the wave function: SSE/Energy. The fraction is known as the Fractional Error.
A theoretical constraint on the above described Gerchberg-Saxton process is that the sum squared error (SSE), and hence the Fractional Error, must decrease or at worst remain constant with each iteration of the process.
Although the Gerchberg-Saxton solution has been widely used in many different contexts, a major problem has been that the algorithm can xe2x80x9clockxe2x80x9d rather than decrease to a sum square error (SSE) of zero. That is to say, the error could remain constant and the wave function, which normally develops with each iteration, would cease to change. The fact that the SSE cannot increase may in this way trap the algorithm""s progress in an xe2x80x9cerror well.xe2x80x9d See Gerchberg, R. xe2x80x9cThe Lock Problem in the Gerchberg Saxton Algorithm for Phase Retrieval,xe2x80x9d Optik, 74, 91 (1986), and Fienup, J. and Wackerman, C. xe2x80x9cPhase retrieval stagnation problems and solutions,xe2x80x9d J. Opt. Soc. Am.A, 3, 1897 (1986). All of the above-identified publications are hereby incorporated by reference for background. Another problem with the method became apparent in one dimensional pictures where non-unique solutions appeared. Furthermore, the algorithm suffers from slow convergence. To date, there are no alternative satisfactory solutions to these problems with the Gerchberg-Saxton method. Accordingly, there is a need for a system and method that can recover wave front phase information without the drawbacks associated with the prior art.
The present invention is driven by an xe2x80x9cerror reductionxe2x80x9d principle and requires a plurality of samples of the wave front from the object being observed. In one aspect, the invention relies on the fact that the back focal plane of a convergent lens on which the scattered wave from the object impinges contains a wave function, which is directly proportional to the Fourier Transform of the object and is therefore directly proportional to the Fourier Transform of the image plane wave function of the object. In the case where the phase difference from one pixel to any of its neighboring pixels only changes slightly, prior art methods were computationally intensive in trying to distinguish between these slight phase differences. Since the actual back focal plane (BFP) wave transforms to the true image in the Image Plane, by the intervention of the drift space between these two planes (mathematically causing the BFP wave to undergo Fourier Transformation yielding the Image Plane wave), in accordance with the invention a very useful relationship is obtained between the measurements in these two conjugate planes. However, other relationships between the waves in these two planes are achievable by changing the phase and/or amplitude distribution in the BFP. In one aspect of the invention, this can be accomplished by using known but physically different phase filters, in the BFP, whose effects on the BFP phase distribution are known. It is noted that there are other physical methods of effectively changing the phase in the BFP (e.g., the use of defocus). The Image Plane wave resulting from this intervention can be very different from the true object wave, consequently yielding new relationships between intensity measurements in these two conjugate planes. The present invention uses several of these new xe2x80x9csynthesizedxe2x80x9d relationships to drastically reduce computation of the reconstructed wave form, to avoid stagnation in the iterating algorithm, and to avoid certain well known ambiguities in the reconstructed wave function.
In one embodiment of the present invention, a random phase filter is inserted into the Back Focal Plane (BFP) of a convergent lens. This phase filter changes the phase for pixels in the BFP in a known way, thereby changing the resulting image in the Image Plane. The phase distribution of the individual pixels in the BFP can be selected randomly, or according to a desired distribution. In alternative embodiments of the invention, conventional convergent and/or divergent lenses can be used as phase filters.
Using the above filter(s), N different sets of amplitude (intensity) data are obtained from the image plane. That is to say, N different images of the object are created in the image plane. It is noted that in an alternative embodiment of the present invention, wave intensities may be recorded in the BFP as well. Next, each of the N intensity images is processed to obtain a xe2x80x9csyntheticxe2x80x9d wave front using the intensity values measured at the Image Plane and phase values that could be random, or may be selected based on prior knowledge. As a practical matter, any initial phase estimate values will work although, for convenience, initially the phase for each complex pixel can be assumed to be zero. The resulting wave function for each of the N images is then inverse Fourier transformed (using standard fast algorithms), and the known phase shift of each of the corresponding BFP filters is subtracted from each pixel. This is done in turn for each of the N images to obtain N estimates of the wave function at the BFP. The resulting BFP estimates are saved for each of the N images. Then, in accordance with a preferred embodiment these BFP estimates are averaged to obtain a single BFP estimate of the complex BFP wave front.
In an alternative embodiment of the present invention, in which BFP intensity data have been measured along with the N IP images, the amplitude of the BFP wave estimate is changed to the measured amplitude distribution at this point in the iterative process. Then for each of the N IP images, the phase shift of its corresponding filter is added in turn to the single BFP estimate and the N different BFP estimates (differing by the N different phase filter effects) are Fourier transformed to generate N estimates of the wave function at the image plane. Each of the N estimates are then corrected using the actually measured amplitude for that particular image. This correction results in an error value. The above process then is repeated in an iterative fashion until the SSE of all N images is sufficiently small for the purposes of the application. In a typical case, less than 1% of the energy of all N images (i.e., the Fractional Error is less than 1%) can be used.
In another important aspect of the invention, it was discovered that absolute stops placed in the illuminating conjugate plane can also be used with only a slight modification of the processing algorithm to unambiguously recover the phase function of a wave front. The same results can also be accomplished by varying the drift space between the two conjugate planes (one containing the intensity/amplitude of the wave front and the other containing the intensity/amplitude of the Fourier Transform of the wave front). Moreover, lossy phase filters and lossy stops can be used in certain practical applications.
Accordingly, in another aspect of the invention, recovery of the phase information can be accomplished using a set of amplitude filters (hereinafter xe2x80x9cstopsxe2x80x9d). Experiments using stops show that they can be used successfully in addition to or as an alternative to the phase filters discussed above with only small modifications of the processing algorithm. In an important practical application, the use of stops can be applied to build a functioning X-ray microscope.
In another aspect of the invention partial or lossy stops or lossy phase filters, or combination thereof can be used in different embodiments for specific practical applications.
In yet another aspect of the invention, instead of using physical stops or phase filters, the desired set of diffraction images can be created by varying the length of the drift space from the specimen plane. Lens-free microscopes, including an X-ray microscope, can be built in accordance with this embodiment as well.
In particular, in one aspect the invention is a method for recovering phase information of a wave front corresponding to a substantially monochromatic coherent radiation, comprising: (a) irradiating a specimen of material with the substantially monochromatic coherent radiation, the specimen being positioned in a first plane; (b) selectively filtering radiation modulated by the specimen according to N pre-determined filtering patterns corresponding to one or more filters, wherein said one or more filters are positioned substantially at the first plane; (c) for each of the N filtering patterns, capturing spatial intensity values for the selectively filtered modulated radiation at a second plane to produce N corresponding intensity distributions, wherein the second plane is a conjugate diffraction plane with respect to the first plane; (d) processing the N intensity distributions captured in the second plane to provide an estimate of the wave front at the first plane, the step of processing comprising correcting the effect of the corresponding filtering patterns; (e) filtering the provided wave front estimate using the N different filtering patterns to obtain N filtered estimates; (f) processing the filtered estimates to produce N estimated intensity distributions at the second plane; and (g) repeating steps (d), (e) and (f) until an error measure associated with the captured and the estimated intensity distributions in the second plane reaches a predetermined threshold. In a preferred embodiment the filtering patterns are phase filtering patterns, amplitude stops or combination thereof. The number of filtering patterns is typically selected dependent on the desired resolution, or can be optimized using information about the underlying wave front. In a specific embodiment, the radiation is X-ray radiation. Further, the method comprises the step of displaying a wave front using recovered phase information.
In another aspect, the invention is an apparatus, comprising: (a) a source of collimated radiation for irradiating a specimen of material positioned in a specimen plane; (b) a plurality of different stops, each one blocking radiation modulated by the irradiated specimen according to a predetermined blocking pattern; (c) one or more sensors capturing for each of the plurality of stops an indication of the intensity distribution of the modulated radiation in a plane that is a conjugate diffraction plane with respect to the specimen plane; and (d) a processor recovering phase information of the wave front of the modulated radiation from the captured intensity distributions and the predetermined blocking patterns imparted by the plurality of stops. In a preferred embodiment the source of radiation is an X-ray source, and the device can be used as an X-ray microscope.
In another aspect, the invention is an apparatus for processing radiation, comprising: (a) a source of collimated radiation for irradiating a specimen positioned in a specimen plane; (b) one or more sensors capturing an indication of the intensity distribution of radiation modulated by the specimen in a plane that is a conjugate diffraction plane with respect to the specimen plane; (c) a motion mechanism changing the distance between the specimen plane and the conjugate diffraction plane, such as to introduce a predetermined phase shift in the modulated radiation; and (d) a processor for recovering phase information of the wave front of the modulated radiation from a plurality of captured intensity distributions obtained using a plurality of predetermined phase shifts introduced by the motion mechanism.
In yet another aspect, the invention is a method for processing substantially monochromatic coherent radiation modulated in a first plane, comprising: (a) capturing N intensity distributions corresponding to the modulated radiation at a second plane, the second plane being conjugate diffraction plane with respect to the first plane, where the captured intensity distributions are obtained by filtering the modulated radiation at the first plane using N different filtering patterns; (b) processing the N intensity distributions captured in the second plane to provide an estimate of the radiation wave front at the first plane, the step of processing comprising correcting the effect of the corresponding filtering patterns; (c) processing the provided estimate of the radiation wave front at the first plane using the N different filtering patterns to compute N estimated intensity distributions in the second plane; (d) computing an error measure corresponding to differences between the captured and the estimated intensity distributions in the second plane; and (e) iteratively repeating steps (b), (c) and (d) until the error measure computed in step (d) drops below a pre-determined threshold.